clc;
% time constant (ms)
tau = 100;
% define sampling frequency (Hz)
fs = 8E3;
fs=100;%处理心电数据


% define signal duration (s)
duration = 0.01;

% get signal length (samples)
N = round( fs*duration );
N=1868


% generate signal time vector (s)
time = [0:N-1]/fs;
TIME = [0:N-1]/fs;

% generate input signals (sinusoids in noise)
data=importdata('/home/appleyuchi/桌面/Matlab/指数平滑/一次指数平滑/result.txt');
M=data;
X=M;

size(X)%注意，平滑前输入的是列向量


% apply exponential smoothing 
Y = expsmooth( X, fs, tau );
Y=Y(2:1868);%修正指数平滑带来的滞后关系
Y(1868)=0;%因为上面只修正了1867个数据，所以最后一个数据随便幅值，对整体波形并不影响

%-----------------------------下面进行峰值恢复-----------------------------------------------

%基线漂移滤波后的峰值所在时刻点的数据
peak_data=importdata('/home/appleyuchi/桌面/Matlab/指数平滑/一次指数平滑/locs.txt');
peak_data=peak_data';%peak_data是心拍所在的时刻位置,并不是心拍数值

final_data=Y
for i=1:length(peak_data)
    sprintf('---------------------开始寻找下标记----------------------------------------')
    peak_data(i);
    
    index=find(abs(TIME-peak_data(i))<1e-5)%寻找时刻中农等于峰值时刻点的下标
    sprintf('---------------------一轮for循环结束，输出峰值数据----------------------------------------')
    final_data(index)%根据下标记获取"心拍数值"
    
    
    %对心拍附近数据进行修复.我们以峰值(心拍)位置为中心，修复周边的修复五个点，即 峰值(心拍)±0.02s范围修复
    final_data(index-2)=final_data(index-2)+20;
    final_data(index-1)=final_data(index-1)+35;
    final_data(index)=final_data(index)+50;
    final_data(index+1)=final_data(index+1)+35;
    final_data(index+2)=final_data(index+2)+20;
    
end    


Y=final_data

%-----------------------------上面进行峰值恢复-----------------------------------------------



% plot the results
figure('Position', [10 10 500 350], 'PaperPositionMode', 'auto', 'Visible', 'on', 'color', 'w'); 

subplot(211); 
plot( time, X ); 
xlabel( 'Time (s)' );
ylabel( 'Amplitude' );
title( 'Noisy signals' );
set( gca, 'box', 'off' );
axis( [min(time) max(time) min(X(:))-0.1*max(X(:)) 1.1*max(X(:))] );

subplot(212); 
plot( time, Y );
xlabel( 'Time (s)' );
ylabel( 'Amplitude' );
title( 'Smoothed signals' );
set( gca, 'box', 'off' );
axis( [min(time) max(time) min(Y(:))-0.1*max(Y(:)) 1.1*max(Y(:))] );

% print results to png file
print( '-dpng', 'expsmooth.png' );
